A global equation of state of two-dimensional hard sphere systems 
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Hard sphere systems in two dimensions are examined for arbitrary density. Simulation results 
are compared to the theoretical predictions for both the low and the high density limit, where the 
system is either disordered or ordered, respectively. The pressure in the system increases with the 
density, except for an intermediate range of volume fractions 0.65 < v < 0.75, where a disorder-order 
phase transition occurs. The proposed global equation of state (which describes the pressure for all 
densities) is applied to the situation of an extremely dense hard sphere gas in a gravitational field 
and shows reasonable agreement with both experimental and numerical data. 
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A hard sphere (HS) system is a simple and tractable 
model for various physical phenomena. It was used to 
examine disorder-order transitions, the glass transition, 
or simple gases and liquids [^|| both theoretically and 
numerically. The theory that describes the behavior of 
rather dilute hard sphere systems is the kinetic theory 
where particles are assumed to be rigid and col- 
lisions take place in zero time. An extension to Boltz- 
mann's low density theory was introduced by Enskog 
1^,^, taking into account excluded volume effects and 
also momentum transport via collisions. In the limit of 
high densities, the cage-effect, where particles are cap- 
tured by their neighbors becomes important and a 
free- volume theory can be formulated Q . There exists 
no theory, to our knowledge, which is valid for the in- 
termediate densities where the system changes from the 
disordered to the ordered state, however, various theo- 
retical approaches were proposed in the last decades, see 
Refs. [PtPHlO|| and references therein. 

When dissipation is added to the HS model, one has 
the simplest version of a granular gas, i.e. the inelastic 
hard sphere (IHS) model. Granular media represent the 
more general class of dissipativc, non-equilibrium, multi- 
particle systems Attempts to describe granu- 
lar media by means of the kinetic theory are usually 
restricted to certain limits like small densities or weak 
dissipation |l^, Q . Also in the case of granular media, 
one has to apply higher order corrections to successfully 
describe the system under more general conditions ||l5| - 
[isll or for multi-particle contacts (l9|. The result of a 
kinetic theory approach is, in the simplest set of 
balance equations for mass, momentum, and energy with 
constitutive expressions for the transport coefficients, de- 
scribing stress, viscosity, heat conduction, and energy dis- 
sipation. In this letter, we focus on the pressure p, the 
isotropic part of the stress in a hard sphere gas - or in a 
granular gas in its elastic limit r — > 1, with the coefficient 
of restitution r. The model system is periodic and two- 
dimensional (2D) with volume V = Ixly, where Ix and ly 
are horizontal and vertical size, respectively. It contains 



N particles with radii a, and masses m which are located 
at positions with velocities Vi. The fraction of the area 
which is covered by particles is denoted as volume frac- 
tion = Nna'^/V. The kinetic energy is £' = J2iLi '"h 
the temperature is defined as T = E/N in 2D and the 
energy density is E/V. For low and intermediate densi- 
ties v < i^c (with the "crystallisation" density Vc at which 
order becomes important) , the kinetic theory leads to an 
expression for the equation of state, i.e. the dimensionless 
excess pressure due to particle interactions 



P := pV/E - 1 = 2i^g{iy) 



(1) 



For an ideal gas with non-interacting particles, one has 
pV/E = 1 so that P = 0; for non-zero densities one 
has P > since the collisions contribute to the momen- 
tum transport and thus to the pressure, the viscosity, 
the heat-conductivity and the dissipation rate. In cases 
with r < 1 the factor 2 can be replaced by 1 -I- r. The 
pair correlation function at contact g(v) accounts for the 
probability that a collision occurs. Typically, 17(1^) is de- 
termined via a virial expansion around low densities and 
one can use 
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where the subscript 4 indicates that the second term is of 
order 1/(1 — i/)^. The first term in Eq. ^ is the simpler, 
frequently used version (72 (i^) introduced by Henderson 
[|l3|,|2^,|l). Note that the expression in Eq. (|) is slightly 
different from the form in Refs. |lj, |l|. The value of 
54, taken at contact, accounts for the excluded volume 
effect and the increase of the collision rate with density. 
At densities larger than ~ 0.7, an ordered triangular 
structure is evidenced p^-|2^. 

One of the unsolved problems concerning an applica- 
tion of the balance equations to a specific boundary value 
problem is the limited range of validity of Eq. (||). Un- 
der realistic conditions with r < 1, the volume fraction 
can take values ly > i^c so that a solution based on 
Eqs. (□) and (g) can be correct up to Vc only. This is 
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even worse, since the "virial" I'gi'y) also occurs in all the 
other transport coefficients. For the same reason, fortu- 
nately, a generalization of the pressure P to all densities 
will thus enter all the other transport coefficients. This is 
why we examine the equation of state at all densities and 
propose a global equation of state which is then tested by 
a comparison with simulation and experiment. 

For the numerical modeling of the system an event 
driven (ED) method |^,|2^ is used. A change in ve- 
locity can occur only at a collision when the standard 
interaction model, based on momentum and energy con- 
servation is used The post-collisional velocities v' 
of two collision partners, in their center of mass reference 
frame, are given in terms of the pre-coUisional velocities 
V, by 2 — ^^1.2 T [{vi — V2) ■ n] n, with the unit vector 
n, pointing along the line connecting the centers of the 
colliding particles. This model can also be extended to 
the more general case of dissipative particles with rough 
surfaces 

The stress tensor inside a test-volume V (whose 
isotropic part is the pressure p) has two contributions, 
one from the convectional transport of mass and thus 
momentum and the other due to collisions and the re- 
lated momentum transport, for details see Refs. p3|,|2^] 
and references therein. The mean pressure is obtained 
from simulations with different volume fractions v in the 
following. 

The equation of state in the dense, ordered phase 
has been calculated by means of a free volume the- 
ory |^,|2^,^, that leads in 2D to the reduced pressure 
Pfv = co/(^'max — t^) — 1, with cq « 1.8137 as obtained 
from our numerical data. Based on the simulation results 
we propose the corrected high density pressure 



and still within about 3 per-cent in the interval 0.65 < 
1^ < 0.75 M. 
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(3) 



where h-i{x) is a fit-polynomial [l -I- cix + c^x^^^ of order 
three, with ci = -0.04 and = 3.25 

What remains to be done is to merge the low density 
pressure P4 and the high density expression (||). To our 
knowledge, no theory exists, which connects these two 
limiting regimes. For various approaches concerning the 
melting and freezing transition see Refs. [^,|9|, p2|^,|29|, 
|30|. Therefore, we propose the global equation of state 



Q = P4 + m(i/)[Pdc„so - P4 
with an empirical merging function 

1 



1 + exp [-(j/ - iyc)/ma] 



(4) 



(5) 



which selects P4 for v ^ i^^ and Pdonsc for j/ 3> i^c: with 
the width of the transition niQ. The parameters i>c = 
0.7006 and mp = 0.0111 lead to qualitative agreement 
between Q and the simulation results - by far better than 
1 per-cent for the purely ordered and disordered regimes. 




0.2 0.4 0.6 



0.8 



FIG. 1. Global, dimensionless equation of state Q, plotted 
against the volume fraction 1/ with logarithmic vertical axis. 
The dashed and dotted lines correspond to P4 and Pdcnsc, 
respectively, see Eqs. (0), and (^. In the inset, simulation 
data (A'^ — 1628, r — 1) are compared with Q. 

When plotting P against the volume fraction 1/ with 
a logarithmic vertical axis in Fig. |l|, the results for the 
different simulations can not be distinguished from the 
theoretical prediction P4 for ly ^ 0.65. For larger volume 
fractions one obtains crystallisation around Vc ~ 0.70 
and the data clearly deviate from P4. In the transition 
regime, the coexistence of fluid and solid phases can be 
obtained. The pressure is strongly reduced due to the 
enhanced free volume in the ordered phase. The reduced 
pressure data eventually diverge at the maximum packing 
fraction in 2D t-™™" = 7r/(2^/3). Note that one has to 
choose the system size such that a triangular lattice fits 
perfectly in the system, i.e. ly/lx = VSh/2w with integer 
w and even h ~ otherwise the maximum possible volume 
fraction can be smaller. Since our simulations are already 
set up on a perfect triangular lattice, the maximum den- 
sity is approached for v — s- fmax- If, instead, the volume 
fraction is increased by increasing the particle size |^,|2^] 
and the simulations are started from random, low den- 
sity configurations, i^max is not reached due to defects in 
the 2D crystal ||2^, ^ . Thus our global equation of state 
represents the high density, small compression-rate limit. 

The transition by itself is also interesting, since we ob- 
tain a hysteresis loop when the density is increased and 
decreased with a finite rate [Q,^, Especially in the 
transition region, the relaxation time is very large, and 
the inflection in the data (see the inset in Fig. P can 
be due to either the finite relaxation time, the finite- 
ness of the system or the initial and boundary conditions 
1^,^. Note that the analytical expression Q allows 
for a straightforward numerical integration of the den- 
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sity profile (see below), since the fit-parameters are cho- 
sen such that the slope of Q is always positive. This is 
a compromise between the quality of the fit on the one 
hand and the numerical treatability of the function on 
the other hand - instabilities are avoided but also mem- 
ory effects are disregarded. 

For a HS system in a gravitational field with the accel- 
eration g in the negative vertical direction, both density- 
and pressure gradient have to be taken into account. 
In the following, we compute analytically the density 
profile for an ideal gas {v < 0.65); the profile for the 
extremely dense gas is computed numerically using the 
global equation of state and is found to be in excellent 
agreement with the numerical ED simulations, where a 
horizontal wall at z = is introduced in a periodic, two- 
dimensional system of width L = lx/{2a) and infinite 
height. The number density n = n{z) ~ N/V is re- 
lated to the volume fraction by n = v{z) / {tto^). Here, 
we briefly sketch how to obtain an analytical solution for 
the density profile, valid at least for low and intermediate 
densities 

Given the equilibrium of forces, the force —Ldp due to 
the pressure gradient at height z compensates the weight 
nmgLdz of the particles in a layer with height dz, so that 
the differential equation dp/ dz = ^nmg has to be solved. 
In the simplest case, the equation of state of an ideal di- 
lute gasp — nT, separation of variables, and the assump- 
tion of a constant temperature, leads to an exponentially 
decreasing density profile ^{z) — i^^exp (— (z — zq)/zt), 
with 1/ < I'd i^iza) and zt — T/{'mg). In a closed 
system, the particle number N is conserved so that inte- 
gration of n over z determines the volume fraction at the 
bottom Vd — Nira^ / (ztL), in the dilute limit. 

In denser situations (0 < v < 0.65) the pressure can be 
expressed as p = nT{l + 2i>g2{i')) [we do not use (74(2^) in 
order to keep the analysis simple], and integration leads 
to an implicit definition of ^{z): 
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2[.g2(i^o)-.92H] (6) 



with the unknown volume fraction vq at zqj again deter- 
mined by the integral over the density. This leads to a 
third order polynomial for vq, which can be solved ana- 
lytically and has at least one real solution. When 
the theoretical density profile in Eq. (|^) is compared with 
numerical simulations, one obtains perfect agreement for 
v < 0.65 [ p3| . Since the functions g2{^') and g4^{v) are 
wrong at larger densities z^, one cannot expect that the 
pressure and the density profile are correct. 

Using the gobal equation of state, Q, from Eq. (Q) 
instead of 2vg2{v), one has to integrate the differential 
equation dp/dz — —nmg numerically withp = nT{l + Q) 
under the constraint that the particle number is a con- 
stant. In Fig. ^, the volume fraction v is plotted against 
the rescaled height z/zt for both theory and simula- 
tions. Simulation parameters are N = 1000, L = 10, and 



ZT/(2a) = 5.85 (open circles) or N = 3000, L = 50, and 
ZT/{2a) = 0.508 (open squares). In addition, we present 
experimental results from vibrated two-dimensional ar- 
rays of small spheres [Q (solid dots), neglecting the fact 
that this situation is weakly dissipative. Both the quali- 
tative and the quantitative behavior of the density pro- 
file is well reproduced by the numerical solution using the 
global equation of state. All solutions belong to one mas- 
ter curve and can be rescaled by a horizontal shift. The 
averaging result is somewhat dependent on the averaging 
procedure ~ we evidence strong coarse-graining effects in 
the dense, ordered regime with densities v > 0.70. Us- 
ing two methods, one tailord for the ordered regime and 
the other for the disordered regime, however, we obtain 
consistent results. 
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FIG. 2. Volume fraction of the hard sphere g func- 
tion of the rescaled height z/zt for different i/o values as given 
in the inset. Lines are the theoretical predictions with increas- 
ing bottom density, from left to right, open symbols are two 
simulations and the solid dots are the experimental data from 
Ref. 0. 

In summary, we tested existing predictions for the 
equation of state of a 2D hard sphere gas of arbitrary 
density via comparison with numerical simulations and 
experimental data. In the dilute case, the particle corre- 
lation at contact and the collision frequency (and thus the 
equation of state) are nicely predicted by the kinetic the- 
ory expressions up to intermediate densities v ~ 0.65. In 
the dense case, the free volume theory for 2D systems can 
be applied to systems with densities larger than v fa 0.75. 
Finally, a merging function is proposed, which connects 
the low and high density regimes, resulting in a diffcrcn- 
tiable global equation of state for the 2D hard sphere gas 
for arbitrary density. 

The equation of state is used to compute analytically 
and numerically the density profile of an elastic, monodis- 
perse HS gas in a gravitational field. For maximum den- 
sities below Vc, the analytical solution works perfectly 
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well, for higher densities, we used a numerical solver 
(MAPLE). The strange shape of the density profile as 
obtained from simulations is nicely reproduced by our 
theory based on the global equation of state, including 
a wiggle a.t V ^ Vc- We remark that the ED simulation 
method parallels the Monte Carlo (MC) method con- 
cerning the particle-particle interactions, but in contrast 
to MC allows for a definition of time and thus for the 
examination of the dynamics. 

The presented results are obtained from homogeneous, 
elastic systems of arbitrary density. The range of ap- 
plicability, however, is much wider. Since already weak 
dissipation can lead to strong inhomogeneities in den- 
sity, temperature, and pressure, the global equation of 
state is a necessary tool to treat effects like clustering, 
surface- waves, pattern formation, or phase-transition and 
-coexistence by means of a continuum theory. In a freely 
cooling "granular gas", for example, clustering leads to 
all densities between « and v « i^max [ p5| . 

The proposed global equation of state is based on a lim- 
ited amount of data from ED simulations. First checks, 
whether our global equation of state still makes sense for 
different particle-size distribution functions are promis- 
ing, however, the crystallisation effect vanishes for strong 
enough polydispersity [ p3| . What remains to be done is 
to find similar global expressions for other transport coef- 
ficients like the viscosity and the heat-conductivity and, 
furthermore, to extend the presented approach to three 
dimensional systems. 
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